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METHOD FOR IMAGING OF PRE-STACK SEISMIC DATA 
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BACKGROUND OF THE INVENTION 

1 . Field of the Invention 

[0001] This invention relates generally to the field of geophysical prospecting. More 
particularly, the invention relates to the field of seismic data processing. Specifically, the 
invention is a method for imaging prestack seismic data for seismic migration. 

2. Description of the Related Art 

[0002] Seismic migration is the process of constructing the reflector surfaces defining the 
subterranean earth layers of interest from the recorded seismic data. Thus, the process of 
migration provides an image of the earth in depth or time. It is intended to account for both 
positioning (kinematic) and amplitude (dynamic) effects associated with the transmission and 
reflection of seismic energy from seismic sources to seismic receivers. Although vertical 
variations in velocity are the most common, lateral variations are also encountered. These 
velocity variations arise from several causes, including differential compaction of the earth 
layers, uplift on the flanks of salt diapers, and variation in depositional dynamics between 
proximal (shaly) and distal (sandy to carbonaceous) shelf locations. 

[0003] It has been recognized for many years in the geophysical processing industry that 
seismic migration should be performed pre-stack and in the depth domain, rather than post-stack, 
in order to obtain optimal, stacked images in areas with structural complexity. In addition, pre- 
stack depth migration offers the advantage of optimally preparing the data for subsequent AVO 
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(Amplitude Versus Offset) analysis. Pre-stack depth migration has traditionally been performed 
through the application of Kirchhoff methods. However, because of recent advances in 
computing hardware and improvements in the design of efficient extrapolators, methods that are 
based on solutions of the one-way wave equation are starting to be used. 

[0004] It has been well established in the literature that wave equation-based methods offer 
the kinematics advantage of implicitly including multipathing, and can thus more accurately 
delineate structures underlying a complex overburden. However, there has been considerably 
less discussion of the dynamic advantages that might (or should) be realized through wave 
equation-based methods. This is not surprising, since Kirchhoff pre-stack depth migration has 
undergone a much longer evolution than its wave equation-based counterparts. Part of this 
evolution has been the development of various factors or strategies to account for geometrical 
divergence and spatial irregularities in acquisition. 

[0005] Poststack migration is often inadequate for areas that are geologically complex, 
because the migration is performed at a late stage, after the seismic data have been affected by 
the NMO (normal moveout) or DMO (dip moveout) corrections and CDP (common depth point) 
stacking. One prestack option is depth migration of common shot gathers, or shot records. This 
method is thus called "shot record" or "shot profile" migration. Shot record migration can give 
more accurate imagining, better preserve dip, and provide accurate prestack amplitude 
information. These properties have made shot record migration one of the more popular 
methods of wave theory migration. 

[0006] Shot record migration is described in Reshef, Moshe and Kosloff, Dan, "Migration of 
common shot gathers", Geophysics, Vol. 51, No. 2 (February, 1986), pp. 324-331. Reshef and 
Kosloff (1986) describe three methods for migrating common-shot records. Discussion of the 
accuracy and efficiency of shot record migration compared to full pre-stack migration is 
contained in the two articles: Jean-Paul Jeannot, 1988, "Full prestack versus shot record 
migration: practical aspects", SEG Expanded Abstracts, pp. 966-968 and C.P.A. Wapenaar and 
A.J. Berkhout, 1987, "Full prestack versus shot record migration", SEG Expanded Abstracts. 
[0007] Seismic migration comprises two steps: (1) wave extrapolation and (2) imaging. The 
present application deals with the second step of imaging in seismic migration. A corresponding 
method of wave extrapolation in seismic migration is described in co-pending U.S. patent 
application, "Method for Downward Extrapolation of Prestack Seismic Data", by the same 

2 

1 



PGS-03-20US 



inventor and assigned to the same assignee as the present patent application. 
[0008] The key elements necessary for an accurate shot record migration may be derived by 
examining the formal representation for the reflection response from a single, flat interface. Shot 
record migration is then seen to be the set of operations necessary to invert for the plane-wave 
reflectivity. After subsequent generalization to the multi-layered case, these elements include (1) 
definition of the source, (2) extrapolation of the wavefields' phases over depth with an 
amplitude-compensating algorithm that preserves vertical energy flux, and (3) application of an 
imaging principle that is based in some way on the definition of reflectivity. 
[0009J Let the upper half-space contain a point, compressional seismic source with 
coordinates x s = (x Sf y St z^) and spectrum S(co), and let it also contain seismic receivers, such as 
hydrophones, with coordinates x r s (x r ,y r ,z r ). Depth is chosen to be positive downward, and for 
simplicity assume that z r = z s . Finally, let the depth of the interface be z\ A far-field expression 
for the pressure of the upgoing (reflected) wavefield, designated by U, is then given by: 



OO CO 

U(x rl y, t z tt t)*(i/cl) \S(o}) exp(-iat) Jexp[ifc x (x, -*,)] 

; (i) 

Jexp[itt, (y r -y s )]R 



exp(-2/<a£ 0 jz'-zj) 



dk y dk x dco 



where Co and §> = 4( c o) are the velocity and vertical slowness in the upper half-space, co is the 
radial frequency, k x and k y are horizontal wavenumbers, and R is the plane-wave reflectivity at 
the reflector. In an acoustic medium, reflectivity R has the simple form: 
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where p 0 is the density in the upper half-space and p It c/, and £ = ^(c } ) are the density, velocity, 
and vertical slowness, respectively, in the lower half-space. Likewise, the downgoing wavefield, 
designated by D, at the reflector depth z 9 has the form: 
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D{k x 9 k y ,z,G))v (// c 2 0 ) S(o>) expHC*, x s + )] 



exp(-i<ag 0 z) 

^0 



(3) 



[0010] Equation (3) provides a representation that is convenient for theoretical analysis, but 
it is not the most practical form for numerical application. Since extrapolation for shot record 
migration is commonly performed in the (co y x) domain, it is more convenient to represent the 
source in the same domain. Away from the source, the wavefield is given by the inverse 
transformation of Equation (3): 



where x ' = (x \y \z ') specifies the observation point for the wavefield in the subsurface. 
[0011] Equation (4) can be used in at least two different ways. First, the wavefield can be 
specified at the bottom of a homogeneous overburden and then downward continued through the 
underlying, heterogeneous structure. Alternatively, in a land setting, a wavefield can be 
specified at some shallow depth using the local velocity, reverse-extrapolated back to the surface 
and then downward continued through the underlying, heterogeneous medium. 
[0012] Thus, the steps needed to invert the recorded seismic data for the plane-wave 
reflectivity, R, can be identified. First, the source wavefield is modeled to the reflector depth, 
perhaps through Equation (4). Second, the upgoing, recorded wavefield (U) is extrapolated to 
the reflector depth z '. Third, the reflectivity R is computed by 



and then inverse transformed back to the spatial domain. 

[0013] Several imaging principles for shot record migration have been reported in the 
literature. Five common methods are reviewed. Kinematic imaging principles only require the 
upgoing, receiver wavefield £/, while dynamic imaging principles depend on some functional 
relationship between both the upgoing and downgoing, source wavefields U and Z>. 



D(x' 9 a>) = S(a>) 



exp(ia)\x'-x s \/c 0 ) 
c 0 2 |x'-x,| 



(4) 



R(k x ,k y ,z' 9 6>) = 



U(k x ,k y9 z' 9 a>) 
D(k x9 k y ,z',G>) 



(5) 
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[0014] A popular, kinematic approach is given by: 
1 n 

R(x) = - £t/(x,<y y ) cxp[io)jt(x)] , (6) 

where x = (x y y,z) correponds to the image point and t(x) is the one-way time from the source to 
the image point as determined from a ray trace. One can expect it to yield results that are 
partially subject to the limitations of ray theory. In order to account for spreading effects, factors 
must be computed from the ray trace and applied in the same manner as for Kirchhoff migration. 
For a further description of this approach, see, for example, the two articles Reshef, Moshe and 
Kosloff, Dan, "Migration of common shot gathers", Geophysics, Vol. 51, No. 2 (February, 
1986), pp. 324-331, or Reshef, Moshe, "Prestack depth imaging of three-dimensional shot 
gathers", Geophysics, Vol. 56, No. 8, (August, 1991), pp. 1158-1163. 
[0015] An average correlation in the (x 9 oo) domain is given by: 

R(x) = -t d [U(x 9 a>j)D\x 9 a> J )] 9 (7) 
n ~[ 

where the superscript "*" denotes complex conjugation. This approach will yield poorly- 
resolved (low frequency), but robust images. However, it does not yield a valid measure of 
reflectivity. An accurate recovery of divergence effects requires that the reflectivity be 
computed as some form of a ratio between the source and receiver wavefields. 
[0016] The classical definition of reflectivity in the (x,co) domain is given by a 
deconvolution: 



1 « U(x 9 a>g) 



This imaging condition is very close to the theoretically correct method, which entails a 
computation in the wavenumber, rather than the spatial, domain. This method can yield 
"floating point errors" where the source illumination is very weak. 
[0017] A normalized deconvolution approach in the (x,co) domain is given by: 
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\D(x,G)j)\ 

where \D\ = (DD*) 1/2 . This algorithm is less sensitive to noise than the deconvolution approach 
in Equation (8). 

[0018] After introducing pre- whitening, Equation (9) becomes: 



h£ |£>(x,^.)| 2 +* 1 ' 



Equation (10) is probably the most commonly used imaging method in the industry. The pre- 
whitening term a stabilizes the solution in the presence of noise. However, it has the undesirable 
side effect of inhibiting an accurate recovery of geometrical divergence. The optimal value for s 
is highly dependent on the data and varies with each shot record and also with spatial position. 
Considerable testing is thus required to find a single, "best compromise" value. It should be 
chosen as small as possible while still yielding a stable, clean-looking solution. For an overview 
of the approaches discussed with reference to Equations (7) - (10), see, for example, Jacobs, 
Allan, 1982, "The prestack migration of profiles", Ph.D. thesis, Stanford Univ., SEP-34, and the 
classic book, Claerbout, Jon, "Imaging the Earth's Interior", Blackwell Scientific Publications, 
Ltd., Oxford, United Kingdom, 1985, 1999. 

[0019] Valenciano, Alejandro, Biondi, Biondo, and Guitton, Antoine, 2002, 
"Multidimensional imaging condition for shot profile migration", Stanford Exploration Project, 
Report SEP-111, June 10, 2002, pp. 71-81, formulate a type of least-squares inverse for the 
imaging condition in shot record migration. A more robust (smoothed) least-squares inverse 
might be formulated by considering data at spatially neighboring locations in computing the 
inverse at location jc. However, their inverses are all performed in the time domain, rather than 
the frequency domain. In the frequency domain, convolution and deconvolution correspond to 
multiplication and division, respectively. In the time domain, deconvolution involves the 
construction of matrices with time-shifted operators. The least-squares inverse formed in 
Valenciano et al. (2002) is based on these matrices. 
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[0020] Valenciano et al. (2002) also propose a simple variation of the pre-whitening factor 
used in Equation (10) above. A binary (0. or 1.) weight can be added to eliminate e when the 
product U D* exceeds some threshold and retain it when U D* is unacceptably small. Results in 
areas of clear imaging would thus not be biased by the pre-whitening, and "floating point errors" 
would be avoided at other locations. 

[0021] Thus, a need exists for a method for imaging prestack seismic data in shot record 
migration that can be fully automated and recover amplitudes in heterogeneous media. 

BRIEF SUMMARY OF THE INVENTION 

[0022] The invention is a method for imaging prestack seismic data. An individual 
reflectivity is calculated for each frequency in the seismic data. Then, a mean reflectivity is 
calculating over the individual reflectivities. A variance is calculated for the set of reflectivities 
versus frequency. A second variance is calculated for the upgoing wavefield, using the mean 
reflectivity. A spatially varying pre-whitening factor is then calculated, using the variance for 
the reflectivities and the variance for the upgoing wavefield. Finally, a single reflectivity is 
calculated at each location, using the spatially varying pre-whitening factor. 

BRIEF DESCRIPTION OF THE DRAWINGS 

[0023] The invention and its advantages may be more easily understood by reference to the 
following detailed description and the attached drawings, in which: 

[0024] FIG. 1 is a flowchart illustrating the processing steps of an embodiment of the method 
of the invention for imaging prestack seismic data; and 

[0025] FIG. 2 is a flowchart illustrating the processing steps of an embodiment of the method 
of the invention for imaging prestack seismic data. 

[0026] While the invention will be described in connection with its preferred embodiments, 
it will be understood that the invention is not limited to these. On the contrary, the invention is 
intended to cover all alternatives, modifications, and equivalents that may be included within the 
scope of the invention, as defined by the appended claims. 
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DETAILED DESCRIPTION OF THE INVENTION 

[0027] The invention is a method for imaging prestack seismic data. The method of the 
invention is illustrated by two embodiments. The first embodiment is an iterative, non-linear 
inversion to solve for an optimum, adaptive, spatially variable pre-whitening factor. The second 
embodiment is a least-squares inversion to determine the optimum reflectivity for the set of all 
frequencies at each depth. 

[0028] At each depth and at each lateral coordinate, a dynamic (amplitude-based) imaging 
principle is applied in order to compute an estimate for the reflection coefficient. One of two 
imaging methods can be used. The first imaging method computes an average over the set of 
results for all "upgoing over downgoing" computations, with an added pre-whitening factor in 
the denominator that automatically adapts to the noise conditions at each location. The 
adaptation is accomplished by iteratively computing variance estimates for the "data" (receiver 
wavefield) and "model" (reflectivity). The second imaging method combines a least-squares 
analysis with an "upgoing over downgoing" computation at each frequency in order to solve for 
a single reflectivity for the entire set of frequencies. 

[0029] The analysis is restricted to the usual assumptions associated with one-way wave 
propagation, i.e., the neglection of multiple reflections and transmission losses. It is also 
assumed that material properties vary smoothly in space. Within the limits of these assumptions, 
shot record migration can adequately recover pre-stack amplitudes in heterogeneous media. 
[0030] Imaging for shot record migration may be developed by examining, for a single 
interface, the exact reflection response in the (x t t) and (k,co) domains. Let x = (x t y,z) be the 
Cartesian coordinates in an infinite, homogeneous, acoustic medium that has velocity c and a 
point compressional source at x s with time dependence S(t). Then the two-way scalar wave 
equation may be written as: 

V 2 P(x,0 - (\/c 2 ) 92p } X9t) = -(4;r/c 2 ) S(x - x 5 ) S(t) . (11) 
dt 2 

The solution to Equation (1 1) is given by: 
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n/ . -S(t-\x-x Ale) 

P(x,i) = ^ (12) 

c x-x. 



where 



|x-xJ = V(*-*i) 2+ (y--V.) a+ (*-*i) 2 - 

[0031] Following the convention that the forward Fourier transforms for a function g(x,y,z,t) 
are given by: 

00 00 00 

g(k x Jc y ,z,o>)= J J jg(x, y, z, t) exp[/ a t - k x x - k y y] dy dx dt , 



the Weyl integral form of the solution in Equation (12) is: 



P(x,a>) = f y incl j S(o>) \cxp[ik x (x-x s )] 



I exp[z k (y - y s )] dk x dk , 

-t 0)4 



where 



Taking an inverse Fourier transform of Equation (13) over a yields: 



(13) 
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P ^ = {/ 4 ^ c 2) fS(«>)exp(-ifl»0 Jexp[i* x (*-x f )] 

00 



(14) 



Equation (14) describes the time-dependent field of the direct arrival that is observed at spatial 
location x, due to the source at x s . 

[0032] Next, a boundary is introduced at depth z ' and the associated reflection response that 
would be recorded by a seismic receiver, such as a hydrophone, is constructed at (x r ,y r ,z r ). For 
the sake of convenience, the depth of the receiver is chosen to be the same as that for the source, 
so z r = z s . From the "reflectivity method", the reflection response can be obtained by modifying 
the kernel in Equation (14) to include the plane-wave reflection coefficient for each downgoing 
wave component. For a single interface at depth z \ the reflectivity function depends only on the 
apparent slownesses p x = k x /co and p y = ky/co, and it is independent of frequency: 



CO CO 



1 1 

£l J 



1 1 

+ — 

V^2 £l> 



(15) 



where & and £ are the vertical slownesses in the upper and lower half spaces, respectively. A 
reflection recording at (x r ,ynZ r =z s ) is thus given by: 



P(x r ,y r ,z s ,t)= V An 2 c i \ J<S(<y) exp(-io)t) \exp[ik x (x r -x s )] 

°f r • / / m n/ k x k y /x exp(-2l'6)^, \z'-Z s \) 

J CO CO 



(16) 



dk x dk y den. 



[0033] Thus, the migration of an individual shot record involves the following three 
operations. First, the source wavefield is forward extrapolated, i.e. modeled, from the depth of 
the source to the imaging depth. Second, the receiver wavefield is reverse extrapolated from the 
surface to the imaging depth. Third, an imaging principle is applied, which may be either 
kinematic or dynamic. Using these three concepts as a guide, the individual wavefields can be 
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identified in the simple, one-layer response. From Equation (13), the forward-extrapolated field 
of the source at the reflector is given by the downgoing wavefield D: 

D(k x ,k yi z\co) = {i/4x 2 c 2 )sW +k y y s )) e *PH<l^ 1) 

The factor "exp[-i(k x x s +k y ys)]" arises from the Fourier transform of the delta-like, lateral position 
of the source. The factor "exp(-/^ 7 |z -z,|)" is the change in phase due to propagation from 
depth z s to depth z\ and the factor "l/fcogj)" is a weight over vertical wavenumber for the 
amplitude of each plane wave. Finally, the factor "S(&)" corresponds to the Fourier transform of 
the source wavelet. 

[0034] Next, using Equation (16), the upgoing wavefield C/that is recorded at depth z r =z s is 
identified as: 

U(k x9 k y ,z s9 a>) = (i/47r 2 c 2 )s(o) exp[-i(k x x,+k y y s )] 

exp(-2/^ | z'-z s \) k x k y (18) 

: a( — , — ,z ). 

cog } co co 

In order to redatum U(k Xt k y ,z s ; co) to the depth of the reflector, U must be extrapolated with 
reversed sign from depth z s to depth z ': 

U(k x 9 k y9 z' 9 a)) = U(k x ,k y9 z s ,co) exp(i fli£ | z' - z 5 |) 

= {i/47r 2 c 2 )s(a>) exp[-i(k x x s + k y y 5 )] (19) 

exp(-ifl>£ \z'-z 5 |) k x k y 

K{ — , — ,z )• 

CO^ CO CO 

[0035] Comparing Equations (17) and (19), the downgoing and upgoing wavefields at the 
reflector thus have the following relationship: 
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U(k x ,k y ,z\co) = D(k x ,k y ,z\6>)R(^,^,z'), (20) 

CO (O 



which yields the imaging condition: 



co co D{k x9 k y9 z 9 co) 



[0036] As expected, the plane-wave reflectivity is thus given by the ratio of the upgoing and 
downgoing wave types. However, this relationship strictly holds only in the apparent slowness 
or wavenumber domains. Inverse Fourier transforms over k X9 and k y yield, to within factors of n 9 

U(x 9 y 9 z\ co) = D(x 9 y 9 z\m)* R( x, y 9 z\co), (22) 

where "*" denotes convolution between U and R over both coordinates x and y. Frequency 
dependence has been retained in Equation (22) in order to accommodate the redundancy afforded 
by measurements at many frequencies. It is common practice to begin with the asymptotic 
relationship: 

U (x, y, z\co) = D(x 9 y, z\co)> R( x 9 y 9 z\co) , (23) 
which yields the alternative imaging condition: 



*( W >) = to^> (24) 
D(x 9 y 9 z' 9 a>) 



[0037] Thus, for homogeneous media, an amplitude-recovering shot record migration can be 
accomplished with the following algorithm, which must be repeated for each depth z '. Construct 
and extrapolate the source wavefield to each depth, z \ by Equation (17). Then, inverse Fourier 
transform over k x and k y to obtain D(x,y,z \ co). Fourier transform the recorded wavefield 
U(x ri y n z r ; co) over x and y and reverse extrapolate to depth z\ Then inverse Fourier transform 
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over k x and k y to obtain D(x,y,z ';&). Compute R(x,y,z ') by Equation (24), either as an average of 
the ratios over all frequencies, or as an average over power-normalized ratios with pre- 
whitening, or by least squares, as discussed above. 

[0038] Equation (21) suggests that imaging should be performed in the wavenumber domain. 
This is inconvenient, since wavefields extrapolated by finite differencing are normally available 
in the spatial domain. Furthermore, Equation (21) is only valid for flat layers, since it implicitly 
includes the assumption of Snell's law for zero dip. 

[0039] The reflection response in the spatial domain can be expanded in terms of wavefronts. 
The first order solution, at high frequency, is proportional to the delta-like behavior of the 
downgoing wavefield. Higher order terms are proportional to Heaviside functions. This analysis 
would show that, to first order, at high frequency, Equation (24) is a good approximation. One 
can still expect some degradation in lateral resolution as a result of this approximation where 
there are rapid lateral changes in reflector properties and, perhaps, for large incidence angles. 
[0040] Thus, asymptotically, one can argue that it is sufficient to compute the ratio in the 
( co,x) domain. Reflectivities can then be computed independently of dip. Since there is a 
redundancy of information provided by several discrete frequencies, a first approximation can be 
to perform an average as follows: 



« t/(x,fi>.)Z>*(x,a>.) 
P\D(x,<Dj)\ 2 + *(x) 



where e(x) is a pre-whitening term, a small number intended to add stability. Equation (25) is 
analogous to the power-normalized classical definition of reflectivity illustrated in Equation (10) 
above, but with a pre-whitening term that can vary spatially. 

[0041] This spatially varying pre-whitening term e(x) is necessary because, for some 
frequencies, the illumination \D\ may lie at the limit of noise. In this case, the unconditioned 
reflectivity ratio is non-physical and may be abnormally large. At some spatial locations, this 
limitation may hold for most or even all of the frequency contributions, which requires that e 
should vary with location. 

[0042] From linear inverse theory, a spatially-variable pre-whitening term (e(x)) can be 
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estimated from an a priori estimate of the variance for the "model" (a R 2 ) and a measured 
variance for the "data" (ay 2 ). Unfortunately, the former estimate is difficult to guess. One can 
get a reasonably good, initial estimate for ct r 2 by first computing individual estimates for R 9 one 
for each frequency, and then computing the mean and standard deviation for this distribution. 
The resulting model variance can then be used, together with ay 2 , to obtain a final estimate for R 
by setting e(x) = au 2 (x)/c R 2 (x) in Equation (10). 

[0043] Thus, the difficulties involved in empirically estimating a spatially varying s(x) can 
be handled by the following automated procedure. FIG. 1 shows a flowchart illustrating the 
processing steps of a first embodiment of the method of the invention for imaging prestack 
seismic data. 

[0044] First, at step 101, a plurality of n frequencies a>j are selected in the prestack seismic 
data. 

[0045] At step 102, a preliminary, constant pre-whitening factor e is selected. The 
preliminary pre-whitening factor does not vary spatially, so it's use in Equation (25) is analogous 
to the simpler estimate of reflectivity in Equation (10). 

[0046] At step 103, a reflectivity R is calculated for each frequency coy selected in step 101 
Preferably, the reflectivities are calculated using Equation (25) with the preliminary, constant 
pre-whitening factor selected in step 102. 

[0047] At step 104, a mean reflectivity <R(x)> is calculated over the distribution of 
reflectivities calculated in step 103. The mean reflectivity is calculated using any appropriate 
methods, as are well known in the art. 

[0048] At step 105, a variance for the reflectivity, <Jr 2 (x), is calculated over the distribution 
of frequencies calculated in step 103. Preferably, the variance for the reflectivity is calculated 
using the following equation: 

[0049] At step 106, a variance for the upgoing wavefield, au 2 (x), is computed using the 
mean reflectivity <R(x)>, computed in step 104. Preferably, the variance for the upgoing 



U(x,a>j) 
D(x,a)j) 



-<*«> 



(26) 
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wavefield is computed using the following equation: 



a* (x) = -L- J (7(1, ® y ) - <*(*)> • Z>(x, )] 2 . (27) 



[0050] At step 107, a spatially varying pre-whitening factor sfo) is calculated using the 
variance for the reflectivity calculated in step 105 and the variance for the upgoing wavefield 
calculated instep 106. Preferably, the pre-whitening factor is calculated by the following 
equation: 

*( X ) = ^7T- (28) 



[0051] At step 108, a reflectivity R(x) is calculated. The reflectivity is calculated by 
reapplying Equation (25) with the spatially varying pre-whitening factor calculated in step 107. 
Thus, the reflectivity is preferably calculated by the following equation: 



m .l± "W-<*«/> (29) 



[0052] At step 109, it is determined if the reflectivity R(x) calculated in step 10 is 
satisfactory. If the answer is no, the reflectivity is not satisfactory, then the process returns to 
step 105 to calculate another reflectivity value. If the answer is yes, the reflectivity is 
satisfactory, then the process ends at step 110. 

[0053] More than one iteration may be necessary, depending on the character of the data. 
Ideally, several iterations could be performed by cycling through Equations (26) - (29). In the 
interest of speed, however, the preferred embodiment is to perform only one iteration. 
[0054] A second embodiment of the method of the invention can also be constructed. The 
imaging principles illustrated in Equations (6) - (10) all represent some average over individual 
estimates of reflectivity R for each frequency cop Each estimate is the slope of a line, with zero 
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intercept, which relates the upgoing to the downgoing wavefield. The average reflectivity is thus 
the average of the slopes. As an alternative, the second embodiment of the method of the 
invention is a least-squares solution for the single slope that minimizes the squared sum of the 
residuals, subject to the constraint that the line has zero intercept. Now, both D and U are treated 
as vectors with respect to frequency. Then the inverse problem for the single parameter R takes 
the form: 

C/(x, Q)j) = D(x, o)j) • R(x) . (30) 
[0055] Equation (30) has the least-squares solution: 



D H (x)-U(x) 

R(X) - D»(x).D(x)' (31) 



where "//" denotes the complex conjugate of the transpose (Hermetian conjugate). Adding a 
pre-whitening factor a, the inverse in Equation (31) becomes: 



RiX) -D»(x).D(x) + £ - (32) 



[0056] Using the vector notation of Equation (30), Equation (32) can be rewritten as: 



m = • (33) 

n j=\ 



In this formulation, it is not necessary to use a spatially- varying pre-whitening factor e(x). 
Equation (33) can be applied either in the space or wavenumber domains. Since this new 
estimate for R is computed as the ratio of sums, rather than as the sum of ratios, one might expect 
fewer problems with conditioning, even in the absence of the pre-whitening factor e. 
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[0057] Thus, for heterogeneous media, an amplitude-recovering shot record migration can be 
accomplished with the following method. FIG. 2 shows a flowchart illustrating the processing 
steps of an embodiment of the method of the invention for imaging prestack seismic data. 
[0058] First in step 201 , a plurality of depths z ' are selected for the prestack seismic data set. 
The depths are selected from any appropriate source known in the art, such as available velocity 
models for the seismic data set. 

[0059] In step 202, a depth z' is selected from the plurality of depths selected in step 201. 
The depths are selected in a systematic fashion, starting from the top depth and proceeding 
downward depth by depth through the seismic data set. 

[0060] In step 203, the source wavefield D(k x ,k y ,z,co) is constructed and extrapolated to the 
depth, z ' selected in step 202'. The resulting source wavefield D(k Xt k y) z \co) may be obtained, for 
example, by Equation (17). 

[0061] In step 204, the recorded wavefield U(x rf y rt z r ; a>) is Fourier transformed over x and y 
and reverse extrapolated to depth z \ The resulting wavefield U(k x ,k y ,z \co) may be obtained, for 
example, by Equation (18) and (19). 

[0062] In step 205, a reflectivity R(x,y,z') is computed at depth z'. The reflectivity may be 
computed, for example, by Equation (21). The reflectivity may be computed either as an average 
of the ratios over all frequencies, as an average over power-normalized ratios with pre- 
whitening, or by least squares. The reflectivity is preferably computed by one of the two 
embodiments of the invention discussed above. The first embodiment of the invention, a power- 
normalized ratio with a spatially varying pre-whitening factor, is discussed with reference to 
Equations (25) - (29) above. The second embodiment of the invention, a least squares method, 
is discussed with reference to Equation (33) above. 

[0063] In step 206, it is determined if any depth z ' remain to be selected in the plurality of 
depths selected in step 201. If the answer is yes, depths z' remain, then the process returns to 
step 202 to select the next depth. If the answer is no, no depths z ' remain, then the process ends 
at step 207. 

[0064] It should be understood that the preceding is merely a detailed description of specific 
embodiments of this invention and that numerous changes, modifications, and alternatives to the 
disclosed embodiments can be made in accordance with the disclosure here without departing 
from the scope of the invention. The preceding description, therefore, is not meant to limit the 
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scope of the invention. Rather, the scope of the invention is to be determined only by the 
appended claims and their equivalents. 
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